Application of wavelets to singular integral scattering equations 
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The use of orthonormal wavelet basis functions for solving singular integral scattering equations 
is investigated. It is shown that these basis functions lead to sparse matrix equations which can 
be solved by iterative techniques. The scaling properties of wavelets are used to derive an efficient 
method for evaluating the singular integrals. The accuracy and efficiency of the wavelet transforms 
is demonstrated by solving the two-body T-matrix equation without partial wave projection. The 
resulting matrix equation which is characteristic of multiparticle integral scattering equations is 
found to provide an efficient method for obtaining accurate approximate solutions to the integral 
equation. These results indicate that wavelet transforms may provide a useful tool for studying 
few-body systems. 
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I. INTRODUCTION 

Few-body systems provide a useful tool for studying 
the dynamics of hadronic systems. The combination of 
short-ranged interactions and finite density means that 
the dynamics of complex hadronic systems can be un- 
derstood by studying the dynamics of few-degree of free- 
dom sub-systems. Few-body systems are simple enough 
to perform nearly complete high-precision measurements 
and to perform ab-initio calculations that are exact to 
within the experimental precision. This clean connec- 
tion between theory and experiment has led to an ex- 
cellent understanding of two-body interactions in low- 
energy nuclear physics, and a good understanding of the 
three-body interactions. 

Our knowledge of low-energy hadronic dynamics is 
largely due to the interplay between experimental and 
computational advances. A complete understanding of 
even the simplest few-hadron system requires measure- 
ments of a complete set of spin observables which have 
small cross sections and require state of the art detectors. 
At the same time, the model calculations with realistic 
interactions are limited by computer speed and memory. 
In addition the equations are either singular or have com- 
plicated boundary conditions which require specialized 
numerical treatments. 

One of the most interesting energy scales is the one 
where the natural choice of few-body degrees of freedom 
changes from nucleons and mesons to sub-nucleon de- 
grees of freedom. The QCD string tension or nucleon size 
suggest that the relevant scale for the onset of this transi- 
tion is about a GeV. A consistent dynamics of hadrons or 
sub-nuclear particles on this scale must be relativistic; a 
Galilean invariant theory cannot simultaneously preserve 
momentum conservation in the lab and center of momen- 
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turn frames if the initial and final reaction products have 
different masses. Relativistic dynamical models are most 
naturally formulated in momentum space. This is due 
to the presence of momentum-dependent Wigner and/or 
Melosh rotations as well as square roots that appear in 
the relationship between energy and momentum. 

Non-relativistic few-body calculations formulated in 
configuration space with local potentials have the ad- 
vantage that the matrices obtained after discretizing the 
dynamical equations are banded, thus reducing the size 
of the numerical calculations. Equivalent momentum- 
space calculations lead to dense matrices of comparable 
dimensions. In addition, the embedding of the two-body 
interactions in the three-body Hilbert space leads to non- 
localities. Realistic relativistic three-body calculations 
are just beginning to be solved [llQl- Numerical meth- 
ods that can reduce the size of these calculations could 
make relativistic calculations of realistic systems more 
tractable. 

In this paper we explore the use of wavelet basis func- 
tions to reduce the size of momentum space scattering 
calculations. The resulting linear system can be can be 
accurately approximated by a linear system with a sparse 
kernel. It is our contention that the use of this sparse 
kernel results in a reduction in the size of the numerical 
calculation that is comparable to the corresponding con- 
figuration space calculations. The advantage is that the 
wavelet methods can be applied in momentum space and 
are not limited to local interactions. 

The long-term goal is to apply wavelet methods to 
solve the relativistic three-body problem. In a previ- 
ous paper l3|, we tested this method to solve the non- 
relativistic Lippmann-Schwinger equation with a Malfliet 
Tjon V potential. In this test problem, the s-wave K- 
matrix was computed. The wavelet method led to a sig- 
nificant reduction in the size of the problem. We found 
that 96% of the matrix elements of the kernel of the in- 
tegral equation could be eliminated leading to an error 
of only a few parts in a million. 



The success of wavelet method in ^ suggests that the 
method should be tested on a more complicated problem. 
In this paper, we test the wavelet method on the same 
problem without using partial waves. This leads to a sin- 
gular two-variable integral equation, which has the same 
number of continuous variables as the three-body Fad- 
deev equations with partial waves. It is simpler than the 
full three-body calculation, but is a much larger calcula- 
tion than was needed in Ref. 3]. In addition, computa- 
tions that employ conventional methods |j] are available 
for comparison. In solving this problem it is necessary 
to address issues involving the storage and computations 
with large matrices. 

One well known use of wavelets is in the data compres- 
sion algorithm used in JPEG files Q- Our motivation 
for applying wavelet methods to scattering problems is 
based on the observation that both a digital photograph 
and a discretized kernel of an integral equation are two- 
dimensional arrays of numbers. If wavelets can reduce 
the size of a digital image, they should have a similar 
effect on the size of the kernel of an integral equation. 

Given the utility of wavelets in digital data processing, 
it is natural to ask why they have not been used ex- 
tensively in numerical computations in scattering. One 
possible reason is because there is a non-trivial learning 
curve that must be overcome for a successful application 
to singular integral equations. A relevant feature is that 
the basis functions have a fractal structure; they are solu- 
tions to a linear renormalization group equation and thus 
have structure on all scales. Numerical techniques that 
exploit the local smoothness of functions do not work ef- 
fectively with functions that have structure on all scales. 

In [3J, we concluded that these limitations could be 
overcome by exploiting the renormalization group trans- 
formation properties of the basis functions in numeri- 
cal computations. These equations were used to com- 
pute moments of the basis functions with polynomials. 
These moments were used to construct efficient quadra- 
ture methods for evaluating overlap integrals. In addi- 
tion, these moments could be combined with the renor- 
malization group equations to perform accurate calcula- 
tions of the types of singular integrals that appear in scat- 
tering problems. A key conclusion of [3| was that wavelet 
methods provide an accurate and effective method for 
solving the scattering equations. In addition, the ex- 
pected reduction in the size of the numerical problem 
could be achieved with minimal loss of precision. 

There are many kinds of wavelets. In J3] we found 
that the Daubechies-3 6] wavelets proved to be the most 
useful for our calculations. Numerical methods based on 
wavelets utilize the existence of two orthogonal bases for 
a model space. The two bases are related by an orthog- 
onal transformation. The first basis, called the father 
function basis, samples the data by averaging on small 
scales. It is the numerical equivalent of a raw digital 
photograph. The orthogonal transformation is generated 
by filtering the coefficients of the father function basis 
into equal numbers of high and low frequency parts. The 



high frequency parts are associated with another type of 
basis function known as the mother function. The same 
filter is again applied only to the to the remaining low 
frequency parts, which are divided into high and low fre- 
quency parts. This is repeated until there is only one 
low frequency coefficient. This orthogonal transforma- 
tion and its inverse can be generated with the same type 
of efficiency as a fast Fourier transform. The new basis 
is called the wavelet basis. 

For the Daubechies-3 wavelets, both sets of basis func- 
tions have compact support. The support of the father 
function basis functions is small and is determined by the 
resolution of the model space. The support of the wavelet 
basis functions is compact, but occurs on all scales be- 
tween the finest resolution and the coarsest resolution. 

The father function for the Daubechies-3 wavelets has 
the property that a finite linear combination of such func- 
tions can locally pointwise represent a polynomial of de- 
gree two or less. Integrals over these polynomials and the 
scaling basis functions can be done exactly and efficiently 
using a one-point quadrature. 

The mother functions have the property that they are 
orthogonal to polynomials of degree two. This means 
that the expansion coefficient for a given mother basis 
function is zero if the function can be well-approximated 
by a polynomial on the support of the basis function. It 
is for this reason that most of the kernel matrix elements 
in this representation are small. Setting these small co- 
efficients to zero is the key approximation that leads to 
sparse matrices. 

Some of the properties that make the Daubechies 
wavelets interesting for numerical computations are 

• The basis functions have compact support. 

• The basis functions are orthonormal. 

• The basis functions can pointwise represent low de- 
gree polynomials 

• The wavelet transform automatically identifies the 
important basis functions. 

• There is a simple one point quadrature rule that is 
exact for low-degree local polynomials. 

• These are accurate methods for computing the sin- 
gular integrals of scattering theory. 

• The basis functions never have to be computed. 

The above list indicates that wavelet bases have many 
advantages in common with spline bases, which have 
proven to be very useful in large few-body calculations. 
Both the spline and wavelet basis functions have com- 
pact support, which allows them to efficiently model local 
structures, both provide pointwise representations of low- 
degree polynomials, both can be easily integrated using 
simple quadrature rules, and both can be accurately in- 
tegrated over the scattering singularity. One feature that 
distinguishes the wavelet method from the spline method 
is that the wavelet transform automatically identifies the 
important basis functions that need to be retained. With 



splines, the regions that have a lot of structure and re- 
quire extra splines need to be identified by hand. This 
is a non-trivial problem in large calculations. The auto- 
matic nature of this step is an important advantage of 
the wavelet method in large calculations. In addition, 
unlike the spline basis functions, the wavelet basis func- 
tions are orthogonal, and the one-point quadrature only 
requires the evaluation of the driving term or kernel at a 
single point to compute matrix elements. This leads to 
numerical approximations that combine the efficiency of 
the collocation method with the stability of the Galerkin 
method. 

In the next section we give an overview of the prop- 
erties of wavelets that are used in our numerical compu- 
tations. Our model problem is defined in section three. 
The methods of section two are used in section four to 
reduce the scattering integral equation in section three 
to an approximate linear system. The transformation to 
a sparse-matrix linear system and the methods used to 
solve the linear equations are discussed in section five. 
The considerations discussed in this section are impor- 
tant for realistic applications. The results of the model 
calculations are discussed and compared to the results 
of partial-wave calculations in section six. Our conclu- 
sions are summarized in section seven. The complex bi- 
conjugate gradient algorithm that was used to solve the 
resulting system of linear equations is outlined in the Ap- 
pendix. 



TABLE I: Scaling Coefficients for Daubechies-3 Wavelets 



ho 


(1 + yiO + \/5 + 2^/10 )/16\/2 


hi 


{5 + VW + 3\/5 + 2yi0 )/16^ 


h2 


(10 - 2^10 + 2\/5 + 2710 )/16^ 


hs 


(10 - 2VT0 - 2\/5 + 2\/T0 )/16\/2 


hi 


{5 + VTO- 3\/5 + 2yi0 )/16^ 


h5 


(1 + \/T0 -\/5 + 2yT0 )/16y2 



where 



gi = {-iyh2K-i-i 



(4) 



The parameter K is the order of the Daubechies wavelet 
and the hi are a unique set of numerical coefficients that 
satisfy certain relations _6] such as orthogonality of basis 
functions. We employ wavelets of order K — 3, hence- 
forth called Daubechies-3 wavelets. The numerical values 
of the hi are given in Table Q] 

Equation Q is the most important in all of wavelet 
analysis, as all the properties of a wavelet basis are de- 
termined by the so-called filter coefficients, hi. A simple 
property that follows from the hi is that the father and 
mother function both have compact support on the in- 
terval (0, 2K—1). All other basis functions are related to 
the primal father and mother by means of dyadic (power 
of two) scale transformations and unit translations. 



II. WAVELET PROPERTIES 

In our work, we use Daubechies' original bases of com- 
pactly supported wavelets [g. In addition to their sim- 
plicity, these functions possess many useful properties for 
numeric calculations, which are discussed at the end of 
this section. 






(5) 



To solve the two-dimensional integral equation for the 
T-matrix we need to construct a two-dimensional direct 
product basis consisting of functions of the form 



A. General Wavelet Analysis 

There are two primal basis functions called the father, 
0, and mother, 1}}. The primal father function is defined 
as the solution of the homogeneous scaling equation 



4'm,l{x)4>n,kiy), 
1pm,lix)(j)n,kiy), 



and 



4>m,li^)'4'n,kiy), 
'lpm,lix)^n,kiy)- (6) 



The primal versions of these four basis function types for 
the Daubechies-3 wavelets are shown in Fig. ^ 



2K-1 



(1) 



(J3{x) = V2 J2 hi(l){2x-l), 
1=0 

with normalization 

{x)dx = 1. (2) 



The primal mother function is defined in terms of the 
father by a similar scaling equation. 



2K-1 



i;{x) = V2 J2 giH'^x-i), 



(3) 



1=0 



B. Equivalent Representations and Wavelet 
Transforms 



If one includes wavelets of all scales, then one can ob- 
tain a basis for L^(M). In practice however, one chooses a 
fine approximation scale J and constructs an approxima- 
tion basis with respect to this scale. At any scale, there 
are two equivalent bases in terms of wavelet functions. 
The first basis consists of translates of the father func- 
tion on the finest scale J. The second basis consists of the 
father functions on the coarsest scale j = and mother 
functions on all intermediate scales j = 0, ..., J — 1. So, 



:^^ ;r^ 



treatment see 0, 13 • First we consider the moments of 
the father function defined by 



{x'') := / x''(f>{x)dx. (9) 

Applying the scahng equation (Q to lO gives 



;^ 




FIG. 1: (Color online) Direct Product Basis of Daubecliies-3 
Wavelets 



for any function we have two equivalent approximations 
given by 



f{x) = ^ai(l)jj{x) 



J-i 



It turns out that the first representation is typically dense 
while the second can often be truncated to a sparse rep- 
resentation by eliminating expansion coefficients with a 
magnitude below some certain threshold. This is because 
the father functions can exactly represent polynomials of 
degree K — 1 while the mother functions arc orthogonal 
to such polynomials 6]. Specifically, 



x'''ip{x)dx = 0, 0<k<K-l 



(8) 



Thus, for any function that is well-represented by low de- 
gree polynomials on the scale J, most of the coefficients 
dj^i in the second representation will be small. These 
small coefficients can be eliminated with a local error of 
0(e), where e is the threshold of the truncation. A fast 
orthogonal transformation known as the discrete wavelet 
transform ^] links the two bases given above. This al- 
lows us to compute projections in the first basis where 
the single scale and single type of basis function make the 
approximations accurate and efficient. Then we can ap- 
ply the discrete wavelet transform to quickly produce the 
sparse basis, which is useful for solving linear systems. 



C. Application of the Scaling Equation 

Now, we briefly discuss some of the useful results that 
follow from the scaling equation Q ■ For a more detailed 



(^ ) ~ ofe 2^ ./o Z^ 



ni—0 



ik—r 



(10) 



This recursion relation, along with the normalization 
condition, (a;") := 1, can be used to compute all of the 
moments of the father function in terms of the filter co- 
efficients, hi. These moments can be used to construct 
quadrature rules, which are used to approximate the pro- 
jection of an arbitrary function, /(x), onto a wavelet 
basis. We employ the simplest such quadrature, the 
one point quadrature |l3. This quadrature is based on 
the identity (x^j = (x) and results in a local error of 
0(/(3)(x)).^ 

It is also important in applications to consider the case 
where the interval of integration is finite. Specifically, 
we consider integrals over left-hand and right-hand end- 
points of the form [Ij 



ix'^y :- 

\ I rn 



[x — m)x dx, 



{x — m)x dx, (11) 



and 



A+ 






{x — m)(j){x — n)dx, 



{x — m)(j){x — n)dx. (12) 



Applying the scaling equation Q to these integrals gives 
linear relations such as 



2K-1 



(-X=2--/2;^ ,,(,.) 



+ 

2m+l 



i=0 



and 



2K-12K-1 



^tn.n = E E ^-^-^ 



+ 
2m+r,2n-\-s ' 



(13) 



(14) 



r=0 s=0 



These linear systems can be solved for the cases of 
m,n — —I, —2, ..., —{2K — 2) using the previously com- 
puted moments for (x*') and the orthogonality relations 
for A+ 

In 13 , we introduced a method for computing singular 
integrals of the form 



Sk-.^ 



4>{x — k) 
x + iO+ 



dx. 



(15) 



TABLE II: Integrals over singularity 



s_l 


-0.1717835441734 


-i 4.041140804162 


S-2 


-1.7516314066967 


+i 1.212142562305 


S-3 


-0.3025942645356 


-J 0.299291822651 


S-4 


-0.3076858066180 


-i 0.013302589081 



where 0"*" is a positive infinitesimal quantity. Applying 
the scaling equation Q), gives the degenerate linear rela- 
tions 



2K-1 



Sk^V2j2 hiS: 



2k-l- 



(16) 



1=0 



These can be supplemented vifith a normalization condi- 
tion 



dx 



/a 
(t){x - n) 
-a 



dx 



x + iO+ 



Y^Sn-.a, (17) 



which follows from the identity 1 — '^^ (p^x — n). Finally, 
we need the nonsingular integrals which can be obtained 
using the recursion relation (|16|l and the convergent ex- 
pansion for large n given by 



^n:a — 



dx 



4>{x — n) 
x + iQ+ 



l + y/n 



dy 



-t(- 



fc^O 



0(j/)y'=dy, (18) 



where the final integrals can be calculated using the 
methods for equations ^ and (|ll|l . The values of the 
singular integrals are given in Table Hll 

For a more thorough and detailed discussion of these 
calculations and additional properties of wavelets see 0, 



III. TWO-BODY T-MATRIX IN MOMENTUM 
SPACE 



The two-body T-matrix is given by the solution to the 
Lippmann-Schwinger equation 



T = V + VGqT, 



(19) 



where V is the two-body potential and Gq ^ {E + ie — 
Hq)^^ is the free two-body propagator. In momentum 
space, this equation becomes 



T{p',p,x') = —vip',p,x\l) - m I dp"p"^ I dx"vip',p",x',x") 



p" 



Po 



:np",p,x"), 



(20) 



where m is the mass of the particles, po is the on-shell momentum, x' = p' ■ p, x" = p" • p, and v is the two- 
body potential with the azimuthal angle dependence integrated out. For our calculations, we use a Malfiiet-Tjon III 
potential |l^ with attractive and repulsive parts. In this case, the azimuthal integration can be carried out analytically 
giving 



v{p',p,x',x) 



1 



Ar 



\/{p'^ + p'^ — 2p'px'x + ^r)^ — 4p'^p^(l — x'^){l — x-^) 

Aa 

\J{p''^ + p^ — 2p'px'x + pa)^ — 4p'^p^(l — a;'^)(l — x'^) 

I 



(21) 



The parameters for this potential arc: Aa — -626.8932 a nuclcon mass such that l/m = 41.47 McV fm^. 
MeV fm, fj.A = 1.55 fm~\ Ar = 1438.723 MeV fm, /iR = In our work, we consider solutions for the half off-shell 

3.11 fm"\ which correspond to those used in Q. We use T-matrix, T{p',po, x'). Traditionally, the T-matrix is de- 



composed in a partial wave basis using 



1=0 



(22) 



where the Pi are Legendre polynomials. Each amplitude 
Ti{p') must be solved for individually. For high ener- 
gies, a significant number of amplitudes may need to be 
included to ensure convergence 4\. 

The magnitude squared of the on-shcU T-matrix is pro- 
portional to the differential cross section. Furthermore, 
the on-shell partial wave amplitudes, Ti{po), can be pa- 
rameterized as 



T/(po) = — — e' 
IT mpQ 



*'(P«)sin(,5i(po)), (23) 



where the Si{po) are experimentally determined phase- 
shifts. These phase-shifts are used to fit realistic nucleon- 
nucleon potentials and should be accurately reproduced 
by any viable solution method. 



IV. WAVELET REPRESENTATION 

To solve equation (|20|l we need to transform the half- 
interval, [0, oo), corresponding to the momentum variable 
into a finite interval, [—a,b]. For computational conve- 
nience we also transform the interval, [—1, 1], associated 
with the angular variable into the region, [— c, d]. For the 
first transformation we use the following map 



b a + k abjp - po) 

p{k):=po-- -, k{p):^ ■ -, 24 

a b — k ap + pob 

which maps the scattering singularity at p" — po to the 
origin. Then we have 



and 



b{b + a) 
'^P-Po-J^^dk 

1 a{b-k) 1 



(25) 



(26) 



p- Po (a + b)po k 
The second mapping is the simple linear transformation 



We now apply these maps to equation (|20() to obtain 
an equivalent integral equation on the rectangular region 
[—a, 6] X [—c,d]. For notational convenience we define 



f{p',x'):=T{p\po,x'), 
g{p',x') := -v{p',po,x',l), 

and for the non-singular part of the kernel 

v{p',p",x',x")p"^ 



L{p',p",x\x") :=TO- 
Now, we let 



p" + Po 



(29) 



(30) 



and 



f{k',u'):^f{p{k'),x{u'j}, 
~g{k',u'):=g{p{k'),x{u')), 



Lik',k",u\u") : 
Lipik'),pik"),x{u'),xiu")) 



(31) 



d + cb-k" 



(32) 



The last factor in this equation comes from applying 
equations (jSSJ, ^^^ and f^ . which gives 



' dpV = l, 2 



P - Po 



k" d + cb- k" 



-dk"du" 



(33) 



Finally, substituting equations H31|) and H32|l into equa- 
tion (l^ni) gives 



fik',u')^g{k',u') 
dk" fdu" ^^^'' ^"]f' ""^ f{k", u") (34) 



— a -' — c 



Now, we project this equation onto the wavelet ba- 
sis which results in a Galerkin type procedure. In gen- 
eral, one can choose a separate fine scale in each variable. 
For notational simplicity, we will consider the case where 
Jk = Ju = J ■ In this case, we approximate / using 



,(,) := ^"-^ + ^ u{x) := ^d + -)-+^d--) ^ (27) 



J{k', U') « Y. fm,n(l>jMk>jAu')- (35) 



which gives 



dx = 



-du. 



Substituting this in 134(1 and multiplying by 
'pj.7n'{k')<pj,n'iu') and integrating over k' and u' 
(28) gives the linear equation 



/ J -^^ m.' .71.' -.m .71. J m .n ~ 9m' ,n' 
m.n 

L{k' ,k" ,u' ,u") 



Y. I dk'f du'l dk"f du"4>j,^,{k')4>j,,,{u') 

mjiJ—a J —c J —a J —c 



4'.Lmik")4>.ln{u")fm,n, 



where 



and 



dk' / du'g{k', u')(f>j^rn' {k')(j>j^n' {u) 



dk' / du'4>j^rn' {k')(f).Ln' («')'/' J,m ik')(t>.],n {u') 



(36) 



(37) 



(38) 



We can evaluate ^m'.rj/ using the one-point quadrature [lOj discussed earher and an endpoint quadrature based on 
the partial moments [2|. Njn',n';m,n is simply the direct product of block diagonal matrices consisting of identity 
blocks and blocks of the form A^ given in equation l|12() . The final term in equation H36(l can be evaluated using the 
subtraction 



-^rn',n';rn,n 



dk' I du i dk" i du"(l)j^rn'{k')(j),j^n'{u) '—^7 '■ 4'.J,m{k")4>.J,n{u") 



'dk'j''du'j\k"j^du"4>j,^, {k')(j>j^^. {U') ^^^'' ^"' "'' ''"\,, ^^^'' "' "^' "''^ ct>J,m{k")cpjAu") 
+ f dk'f du'f du"cj)j,^,{k')4>J^n'{u')L{k',Q,u',u")4>J,n{'^") f '^^^^dk" 



— a J —c J —c 



(39) 



r 



The first term in this equation is nonsingular and can 
be approximated using the quadrature methods previ- 
ously discussed. Likewise, the k',u",u' integrations in 
the second term can be carried out in the same manner. 
The final integration over k" can be accomplished using 
the method following equation p6|) . 

Thus, the problem is reduced to solving a linear system 
of the form 



m.n 



iN„ 



+ L 



m' ,n' :m,n * -^m' ,n' im.nj J m.n 



)U 



(40) 



Once we have solved this equation for fm.n we can substi- 
tute this approximate solution back into the right hand 
side of the original equation to obtain a refined solution. 



V. WAVELET TRANSFORM AND SPARSE 
SOLUTION 

The eigenvalues of A+ accumulate at while those of 
A~ accumulate at 1 as iiT increases [l^. This makes the 
matrix N, and consequently the right hand side of (|40|l . 
numerically ill-behaved. To circumvent this difficulty we 
can precondition the system by inverting N, which is 
easily accomplished by inverting the two blocks A"*" and 



A and using the direct product structure of N. If we 
define 



h = Nf, 

then equation (|40|l becomes 



(I + LN-i)h = g. 



(41) 



(42) 



If we define 



A=(I + LN-i). (43) 

Then equation (|42|l is a simply linear system of the form 



Ah = g. 



(44) 



This is a large dense linear system. However, as shown 
in equation Q, there are two equivalent representa- 
tions that are linked by a fast orthogonal transforma- 
tion. In two variables, the matrix representation of this 
transformation is simply the direct product of the one- 
dimensional transformation matrices that are given in 



many standard references |3- If we denote this matrix 
as W then we can transform equation (|44ll as 



(WAW'^)Wh = Wg. (45) 

Now we make the definitions 

A = WAW'^, h = Wh, g = Wg. (46) 

Then as mentioned in reference to equation Q we can 
truncate the matrix A by ehminating all elements with 
a magnitude below some certain threshold e, where the 
error introduced is proportional to e. The matrix A can 
be stored in a smrse format such as compressed column 
format (CCS) [lj|i which permits both efficient storage 
and matrix multiplication. These savings help eliminate 
computationally costly writing and reading from the hard 
disk when solving the linear system. 

To solve the sparse system we use the complex bicon- 
jugate gradient method 0,113 which we present for gen- 
eral complex matrices in Appendix IXI This method is a 
simple and effective iterative method for general sparse 
matrices. In addition, this method, like all iterative tech- 
niques, is readily amenable to parallel processing of the 
matrix multiplication. Once the solution to H45I) is found, 
it is a simple matter to recover f by applying the inverse 
transform W'^ and the inverse matrix N^^. In particu- 
lar, 



f = N"^W^h. 



VI. RESULTS AND ANALYSIS 



(47) 



We made calculations using the Daubechies-3 wavelets 
at scales up to J = 5 in each variable. The total number 
of wavelet basis functions in each variable was taken to 
be 2*^, where M = J + 2. If we take a = c = 1, then h 
and d are determined by 



h = 2" •^'=(2*=^'= - (2i^ - 2)) - a, 
d = 2--'-{2''^- -{2K-2))-c. 



(48) 



Using these parameters calculations were performed at 
lab energies of 300 and 800 MeV to test the efficacy of the 
method in different energy regimes. Figure [21 shows the 
real and imaginary parts of the half off-shell T-matrix 
as a function of momentum, p', and the scattering an- 
gle, x' = cos{9) at a scattering energy of 800 MeV. 
Daubechies-3 wavelets were used with M^ = M„ — 5. 

It can be seen that the real part of the T-matrix is 
relatively smooth, while the imaginary part does have 
some structure. The fact that the T-matrix is smooth 
with isolated structure suggests that wavelet methods 



Re(T) 



lm(T) 




20 1 




FIG. 2: (Color online) Half off-shell T-matrix at 800 MeV 



should be able to efficiently compress the matrix A. Fig- 
ure|31shows the on-shell T-matrix as a function of angle, 
x' = cos(^), at a scattering energy of 300 MeV. These 
calculations were made using Daubechies-3 wavelets with 
Mk = Mu = 5. From these graphs the general smooth- 



Re(T) 



lm(T) 





FIG. 3: On-shell T-matrix at 300 MeV 



ness of the on-shell amplitude is apparent. We can also 
see the forward peaking of the scattering amplitude that 
is expected at higher energies. 

From a numerical standpoint, the first aspect of the 
calculation to consider is the general convergence of the 
method as the number of basis functions is increased. Ta- 
bles mil and IIVI illustrate the convergence of the method 
as the number of basis functions is increased. The quoted 
values are for on-shell scattering at an angle of 90° using 
Daubechies-3 wavelets with no truncation. From Table 
mil we see that the majority of improvement occurs as 
Mk is increased. This can be attributed to the fact that 
the integral over k" in the kernel is singular and thus 
requires more basis functions to accurately represent the 
dependence on this variable. 

All of these calculations were made by solving the lin- 
ear system (I42f) . It is instructive to consider the behav- 



TABLE III: Convergence as a function of total number of 
basis functions: 300 MeV 



Mk 


Mu 


Re(r(po,po,0)) 


Im(r(po,Po,0)) 


4 


4 


0.484065410 


0.292438234 


4 


5 


0.483906981 


0.293504057 


5 


4 


0.491111143 


0.286418162 


5 


5 


0.490972783 


0.287464852 


5 


6 


0.490891484 


0.287452418 


6 


5 


0.491773044 


0.286276262 


6 


6 


0.491691220 


0.286263888 


6 


7 


0.491678773 


0.286256958 


7 


6 


0.491772680 


0.286123199 


7 


7 


0.491760404 


0.286116271 



TABLE IV: Convergence as a function of total number of 
basis functions: 800 MeV 

'W Mu Re(r(po,po,0)) Im(r(po,po,0)) 



4 


4 


0.456127838 


0.126626540 


4 


5 


0.454750689 


0.126515462 


5 


4 


0.456227507 


0.113862313 


5 


5 


0.455006434 


0.113967425 


5 


6 


0.454697067 


0.113367584 


6 


5 


0.455242066 


0.111684819 


6 


6 


0.454931562 


0.111107815 


6 


7 


0.454884387 


0.111005571 


7 


6 


0.454978565 


0.110889988 


7 


7 


0.454931334 


0.110788571 



ior if one attempts to solve H40() using iterative methods 
without directly inverting N first. Table compares 

the error in the residual, e„ = ||f„|| = g — Ah„ , as a 

function of the number of iterations. As the number of 
iterations, n, is increased the preconditioned method con- 
verges very rapidly, while the non-preconditioned method 
fails to converge adequately. 

Now we turn our attention to the compression of the 
sparse matrix and its subsequent effect on the calculation. 
Figure 01 displays a such a representation for scattering 
at 800 MeV using Daubechies-3 wavelets with Mk = 4, 
Mu = 3. The plot shows the location of the nonzero 
elements of A after it has been truncated at the thresh- 
old level e = 10~^. This threshold produces a matrix 
with 19% of the elements of the full matrix. The order- 
ing scheme for A used in the plot places the elements 
associated with finer scales at higher indices. The de- 
gree of sparsity increases considerably as the scale in- 
creases, which demonstrates that less and less elements 
are needed at finer scales. For the range of test cases we 



TABLE V: Convergence of the bi-conjugate gradient method 
n en (non-preconditioned) e,i (preconditioned) 

1.98xl0~'' 
3.07X10"-""' 
2.03x10"^ 
9.10x10"^" 
8.09x10"" 



10 


3.25x10 


20 


4.21x10 


30 


7.23x10 


40 


9.94x10 


50 


1.01x10 
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FIG. 4: Location of the nonzero of elements of A 



considered, the solution of the sparse linear systems only 
took 5-10% of the computational time. The majority 
of the time was spent calculating and filtering the dense 
matrices. 

In Table IVII the effect of truncating the matrix A on 
the convergence of the solution is illustrated for various 
threshold levels with a lab energy of 800 MeV. This cal- 
culation was performed using the Daubechies-3 wavelets 
with Mk = Mu = 6. Comparing these results with those 
in Table IIIII we see that even keeping just one percent 
of the matrix elements we are able to reproduce the T- 
matrix to the same precision as the accuracy of the un- 
truncated matrix. 

Finally, we consider the accuracy of the phase shifts 
determined by our momentum vector approach. To cal- 
culate the phase shifts we project our T-matrix onto the 
partial waves using 



Ti{p') = 2nJ Pi{x')T{p',po,x')dx'. 



(49) 



We compute the integrals using 20 Gauss-Legendre 
points 16]. From the T/ it is straightforward to calculate 
the phase shifts using equation 123|) . Table IVTll displavs 
these phase shifts (calculated for i?Lab = 800 MeV using 
Daubechies-3 wavelets with Mk = M„ — 7 and truncat- 
ing all but 2% of the coefficients) compared with phase 
shifts calculated using standard partial wave techniques. 
The agreement between the two methods is very good for 
all the phase shifts. 
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TABLE VI: Effect of truncation on the on-slieli T-matrix at 
SOOMeV for scattering at 180°, 90° and 0° corresponding to 
T{po,po,-l),T{po,po,0) andr(po,Po,+l) 



e 


% 


Re(r(po,Po, 


-1)) Re(r(po,po: 


-1)) 





100 


0.249235 


-0.0777091 




io-« 


23 


0.249235 


-0.0777093 




10-'' 


14 


0.249234 


-0.0777116 




io-« 


8 


0.249217 


-0.0777525 




lO"'^ 


1 


0.248296 


-0.0770660 




e 


% 


Re(r(po,Po, 


0)) Re(r(po,PO: 


.0)) 





100 


0.454932 


0.111108 




io-« 


23 


0.454932 


0.111108 




10-'' 


14 


0.454932 


0.111108 




lO"** 


8 


0.454941 


0.111117 




10"^ 


1 


0.454966 


0.111154 




e 


% 


Re(r(po,Po, 


+1)) Re(r(po,PO: 


.+1)) 





100 


-6.16347 


-1.31548 




io-« 


23 


-6.16347 


-1.31548 




10-' 


14 


-6.16347 


-1.31548 




10"" 


8 


-6.16346 


-1.31548 




10"^ 


1 


-6.16327 


-1.31559 





TABLE VII: Comparison of 800 MeV phase shifts with stan- 



dard methods 



I 


Si{po) (Standard) 


Si{po) (Wavelet) 





-0.2535 


-0.2534 


1 


0.2950 


0.2949 


2 


0.3635 


0.3634 


3 


0.2747 


0.2746 


4 


0.1755 


0.1755 


5 


0.1053 


0.1052 


6 


0.06169 


0.06168 


7 


0.03591 


0.03591 


8 


0.02089 


0.02089 


9 


0.01217 


0.01217 


10 


0.007110 


0.007109 


11 


0.004164 


0.004163 


12 


0.002445 


0.002444 


13 


0.001439 


0.001437 



VII. CONCLUSIONS 

We have shown that it is possible to use wavelets 
to calculate the two-body scattering matrix in terms of 
momentum vectors without resorting to partial waves. 
We were able to accurately reproduce the phase shifts 
of the Malfliet-Tjon potential. These calculations lead 
to sparse matrices, which can be efficiently inverted us- 
ing standard iterative methods. Application of a sim- 
ple preconditioning matrix was shown to be necessary 
to achieve convergence of the iterative methods. Tradi- 
tional methods for solving scattering equations in mo- 
mentum space typically produce dense matrices that re- 
quire a large amount of storage and are time consuming 
to invert. These are promising results because relativistic 
scattering equations are naturally formulated in momen- 
tum space. Also, the scattering boundary conditions are 
most easily treated in momentum space. Wavelet meth- 



ods can help treat both of these of problems. 

One of the main advantages of wavelet methods over 
methods such as splines is that the wavelet transform 
presents a method that automatically determines what 
basis functions are necessary for a given accuracy. Unfor- 
tunately, this also leads to one of the main drawbacks of 
this method. In our procedure, a large dense matrix. A, 
needs to be produced first and then this is transformed 
to a sparse matrix. Most of the computational time is 
spent constructing and transforming this matrix into a 
sparse format. The subsequent solution of the sparse lin- 
ear system takes relatively little computational effort. 

For this specific problem, wavelet methods based on 
momentum vectors may not be necessary. The maxi- 
mum number of partial waves that needs to be included 
to achieve convergence, Imax = 14 |j], is simply too small 
to gain a computational benefit from using wavelets in 
the angular variable. To achieve a computational benefit 
we should use less basis functions in the angular vari- 
able than the maximum number of partial waves. In the 
three-body problem or at much higher energies, the num- 
ber of partial waves that need to be included increases 
considerably and computational benefits may be gained 
from employing a momentum vector approach. 
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APPENDIX A: COMPLEX BICONJUGATE 
GRADIENT METHOD 

The biconjugate gradient method [3)ll3 i^ an iterative 
technique for solving large matrix equations of the form 



Ax: 



The advantage of this method for large sparse matrices is 
that it only involves matrix multiplication by A and its 
adjoint, both of which can be accomplished efficiently in 
a sparse storage format such as CCS 14]. The algorithm 
generates a sequence of approximate solutions, x^^ with 
residual r^ = b — Ax^. One iterates until the norm of 
the residual is less than some predetermined value. 

This method is traditionally formulated for real matri- 
ces, but the extension to complex matrices is straightfor- 
ward. Below we present the algorithm for general com- 
plex matrices. For our calculations, we start with the 
initial approximate solution 



with the residual 



xo = b 



ro = b - Axo . 
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For the initial values of the bi- residual Vq, the direction 
vector po, and bi-direction po we use 

ro = b - A'fxo 

Pi =ro 
Pi = Fo . 

Then we use the recurrence relations 



ttfc = 



^Li. 



plAPfc 
Xfe = Xfc-i +akPk 



Tfe 






st 



Pfe+i = r/c + fikPk 
Pfc+i = Ffc + /3fcPfc 



to generate an improved approximation. This is repeated 
until the desired accuracy is obtained. We measure the 
accuracy by the i"^ (C) norm of the residual, 



efc = llrfcll = 
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